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£pH ABSTRACT 

Very-Mgh-energy (VHE; E > 100 GeV) and high-energy (HE; 100 MeV < E < 100 GeV) data from y- 
,H ray observations performed with the H.E.S.S. telescope array and the Fermi-LAT instrument, respectively, are 

analysed in order to investigate the non-thermal processes in the starburst galaxy NGC253. The VHE y-ray 
data can be described by a power law in energy with differential photon index T = 2.14 ± 0.18 stat ± 0.30 sys and 
differential flux normalisation at 1 TeV of F = (9.6 ± 1.5 stat (+5.7, -2.9) sys ) x 10~ 14 TeV _1 cirT 2 s -1 . Apower- 
£NJ law fit to the differential HEy-ray spectrum reveals a photon index of T = 2.24±0.14 stat +0.03 sys and an integral 

flux between 200 MeV and 200 GeV of J F(0.2-200 GeV) = (4.9+ 1.0 stat + 0.3 sys )x 10~ 9 cm" 2 s _1 . No evidence 
for a spectral break or turnover is found over the dynamic range of both the LAT instrument and the H.E.S.S. 
experiment: a combined fit of a power law to the HE and VHE y-ray data results in a differential photon index 
T = 2.34 + 0.03 with a p-value of 30%. The y-ray observations indicate that at least about 20% of the energy 
of the cosmic rays capable of producing hadronic interactions is channeled into pion production. The smooth 
alignment between the spectra in the HE and VHE y-ray domain suggests that the same transport processes 
dominate in the entire energy range. Advection is most likely responsible for charged particle removal from 
CN the starburst nucleus from GeV to multiple TeV energies. In a hadronic scenario for the y-ray production, the 

single overall power-law spectrum observed would therefore correspond to the mean energy spectrum produced 
by the ensemble of cosmic -ray sources in the starburst region. 
• »h Subject headings: Galaxies: starburst, Galaxies: individual: NGC 253, Gamma rays: galaxies, Radiation 

/\ mechanisms: non-thermal, Diffusion, Advection 
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1. INTRODUCTION 

Starburst galaxies are galaxies that undergo an epoch of 
star formation in a very localised region (the starburst re- 
gion) at a rate that is enhanced in comparison to other, so- 
called late-type galaxies such as the Milky Way galaxy. It 
is believed that this starburst activity is triggered either by 
galaxy mergers, a close fly-by of galaxies, or by Galactic 
bar instabilities, where the dynamical equilibrium of the in- 
terstellar gas gets disturbed. This leads to the formation of 
regions of very high-density gas, usually at the centre of the 
galaxy, and subsequently to star formation and a strongly in- 
creased supernova (SN) explosion rate. SN remnant shocks 
are widely believed to be acceleration sites of cosmic rays 
(CRs). This is one reason why starburst regions might have 
a high CR density. Given the high density of target mate- 
rial that is available for p-p interactions and the production of 
7r°s, the starburst nucleus is in addition a promising source of 
high-energy (HE; 100 MeV < E < 100 GeV) and very-high- 
energy (VHE; E > 100 GeV) y rays. From energetic elec- 
trons also Bremsstrahlung and Inverse Compton y rays are 
expected. These electrons may be either directly accelerated 
by the same processes as the nuclear particles, or be gener- 
ated in the decays of charged pions from hadronic collisions. 
They might, however, also be produced in different sources, 
like in pulsar wind nebulae. Starburst galaxies have been pre- 
dict ed early on to be detectable by present y-ray instruments 
(e.g. |Volk et al.|1989|[l996l|Akyuz et al.|1991||Paglione et al. 
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1996). 

The spiral galaxy NGC253 is the closest object in the 
southern sky that belongs to the class of starburst galaxies. 
Its distance is measured as (2 .6-3.9 ) Mpc usi ng different dis- 
tance estimation techniques ([D avidge et al.|1991| |Karachent>] 
|sev et al.]|2003| |Rekola et al.||2005|l. The reference distance 
is d — 2.6 Mpc ( |Davidge et al.|199l) since this value is used 
most widely in the literature to determine the properties of 
NGC253. However, this reference dista nce has recentl y been 
convincingly revised to 3.5 Mpc (Dalcanton et al. 2009). The 
final numerical values used below will therefore be needed to 
be scaled for consistency with the revised distance value. 

Compared to the Milky Way galaxy, NGC 253 exhibits an 
increased overall star formation rate (SFR), with the SFR in 
the starburst nucleus being comparable to that in the entire re- 
maining disk of the galaxy. A SN rate vsn of this nucleus can 
be determined from the far infrared (FIR) observations, since 
the FIR lumin osity can be assumed to be directly propor- 
tional to v S N ( |Van Buren & Greenhouse|1994| >. For NGC 253 
as a whole the SN rate is estimate d to be » 0.08 yr" 1 , with 
» 0.03 yr -1 in the starburst region ( Engelbra chtetafll 998 >■ 
By assessing the SFR |Melo et alT] ( |2002| l found that it can 
amount to 5M yr~' in the starburst nucleus alone which 
is 70% of the SFR of the entire galaxy. The starburst re- 
gion itself has a cylindrical shape with a radius of m 150 pc 
and a full height of =s 60 pc perpendicular to the disk of 
the galaxy and symmetric to its mid-plane with a v olume 
V SB * 1.2 x 10 62 (c//2.6Mpc) 3 cm 3 < |Weaver et al.|2002 1. 

To understand the observed y-ray emission, a simplified 
scenario is considered in which the y-ray production result- 
ing from particle acceleration in the part of the disk outside 
the starburst region is neglected in comparison with that from 
the starburst region. The reasons are the low average gas den- 
sity and radiation field intensity of the average Interstellar 
Medium, and the expected dominance of energy-dependent 
diffusive particle losses from the disk - quantitatively similar 
to the situation in the Milky Way g alaxy. This ex pectation is 
also consistent with the estimate of Strong et al. (2010), who 
find that the HE y-ray luminosity of the Milky Way galaxy is 
an order of magnitude lower than the y-ray luminosity of the 
starburst region of NGC 253 (cf. Sectional. 

The stellar winds from the early-type stars and the sub- 
sequent core collapse SN explosions heat the lower-density 
parts of the surrounding material, causing them to expand 
rapidly from the starburst region in the form of a collective 
wind. The shocks from the SN explosions are the primary 
accelerators of CRs in this scenario and their pressure adds 
to the excess thermal gas pressure. The dense material in the 
starburst region outside the SN remnants will remain essen- 
tially non-ionised in this process and will not participate in 
the flow. The 7r°-producing CRs from the percolating wind 
flow are nevertheless likely to penetrate also the dense gas 
in the starburst region which therefore is a massi ve target for 
y-ray p roduction. The recent detections of HE ( |Abdo et al.| 
2010a) and VHE y-ray emission from the starburst galaxies 
NGT253 dAcero et al .12009) and M 82 ( |Acciari et al.|2009) 
appear to support this picture. 

In a picture where there is quasi-steady equilibrium be- 
tween production and loss processes, the population of high- 
energy CRs accelerated in NGC253's starburst nucleus is re- 
moved from the starburst region predominantly via three dif- 
ferent processes: (i) advective removal of particles in the star- 



burst wind; often also called a "superwind" (see e.g. Weaver 
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[et~aT1|2002l |Zirakashvili &"V 61k 2006), not to be confused 
with the large-scale galactic "disk wind" (see e.g. Breitschw 



erdt et al.[| l991 ; 



Heesen et al.| 2009[ ) that is primarily driven 
by the general population of CRs and the hot gas, both pro- 
duced in the galactic disk, (ii) diffusion of particles from the 
source region, and (iii) catastrophic inelastic ("p-p") interac- 
tions. Energetic electrons/positrons suffer, in addition, ra- 
diative losses. The contributions of these components and 
the resulting y-ray spectra have been discussed by several 



groups (see e.g . Paglione et al. 1996; Aharonian et al. 
Domingo-Santamarfa & Torres 2005 



2005 



2010 



Rephaeli et al. 

Lacki et al.||2010| |2011[ ) and are compared to the measure 
ments presented here. 

After the discove ry of VHE y-ray emission from NGC 253 
(Acero et al. 2009), we present for the first time a spectral 
analysis of the VHE y-ray data obtained by H.E.S.S. in con- 
junction with the analysis of a 30-month set of Fermi-LAT 
data, increased in size by a factor of a; 3 compared to the one 
used in the original publication of the Fermi Collaboration on 
NGC 253 ( |Abdo et al.|2010a| ). These results are used to esti- 
mate the properties of the underlying CR population, such as 
the particle energy density, as well as to place constraints on 
CR transport and the r.m.s. magnetic field strength within the 
starburst region. 

2. H.E.S.S. OBSERVATIONS AND DATA ANALYSIS 

2.1. H.E.S.S. instrument 

The High Energy Stereoscopic System (H.E.S.S.) is an ar- 
ray of four imaging atmospheric Cherenkov telescopes lo- 
cated in the Khomas Highland of Namibia, 1800 m above sea 
level. The telescopes are identical in construction and each 
one comprises a 107 m 2 optical reflector composed of seg- 
mented spherical mirrors and a camera built of 960 photo- 
multiplier tubes. H.E.S.S. ut ilises the im aging atmospheric 
Cherenkov technique (see e.g.; Hillas|1985 i. Cherenkov light, 
emitted by the highly relativistic charged particles in exten- 
sive air showers, is imaged by the mirrors onto the camera. 
A single shower can be recorded by multiple telescopes under 
different viewing angles, allowing stereoscopic reconstruction 
of the primary particle direction and energy with an average 
energy resoluti on of 15% and an ev ent-by-event spatial reso- 
lution of 0.1° (|^aronian[etal]2006}. 

2.2. Data Set 

NGC 253 was observed with the H.E.S.S. array in 2005 and 
from 2007 to 2009 for a total of 241 hours. After standard data 
quality selection, where data taken under unstable weather 
conditions or with malfunctioning hardware have been ex- 
cluded, the total live time amounts to 177 hours of three- and 
four-telescope observations that were used for the generation 
of sky maps of the y-ray emission and the reconstruction of 
energy spectra. Observations were carried out at zenith angles 
of 1° to 42°, with a mean value of 12°. Observations have 
been performed in the wobble-mode, where the telescopes 
were alternately pointed offset in RA and Dec from N GC 253 
resulting in an average pointing offset of 0.5°( Aharonian et al. 
[20061 ). 

2.3. Data Analysis 

All results presented in the following were obtained us- 
ing the Model Analysis (MA; |de Naurois & Rolland||2009 1 
for event reconstruction and background reduction and were 
cross-checked with the boosted-decision-tree-based (BDT) 
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Fig. 1. — Smoothed H.E.S.S. y-ray excess map in units of VHE y-ray events 
per arcmin 2 of the 1.5° X 1.5° FoV, centered on the position of NGC 253. The 
image was smoothed with a Gaussian kernel of 3.9' r.m.s., the radius that 
corresponds to the PSF for this analysis. The black star marks the position of 
the optical centre of NGC 253 and the inlay represents the size of a point-like 
source as it would have been seen by H.E.S.S. for this analysis. White con- 
tours depict the optical emission from the whole galaxy with contour levels 
of consta nt surfac e brightness of 25 mag arcsec~~ and 23.94 mag arcsec~~ 
as used in Pence ( 19801. The dashed circle indicates the 95% error contour 
of the best-fit position of the Fermi-LAT source (see also Table[3j. (Colour 
version of this figure available online.) 



Hillas parameter technique described in detail in Ohm et al. 
(2009). Two different sets of cuts were used for MA: the 
standard cuts require a minimum shower image intensity of 
60 p. e. in each camera. They maximise the acceptance of y- 
ray like events (at the expense of a larger background) and 
are therefore used for energy spectra and sky maps; the faint 
cuts, requiring a higher minimum intensity of 120 p.e., result- 
ing in improved angular resolution at the cost of lower y-ray 
acceptance, are used for position and extension determination. 
Spectral results were derived using the Reflected background 
model, whereas the Ring background model was utilised to 
generate sky maps (Berge Tet al |2007) ). The analysis thresh- 
olds for the standard and faint cuts configuration are 190 GeV 
and 250 GeV for the MA method, respectively. 

2.4. Results 

A VHE y-ray excess map for the 1.5° x 1.5° Field-of-View 
(FoV) centred on the optical position of NGC 253 and pro- 
duced with MA, standard cuts, is shown in Fig. [T] The map 
has been smoothed with a 2D Gaussian kernel of3.9' r.m.s. 
to reduce the effect of statistical fluctuations and matched to 
the PSF for this analysis. A total of 329 + 49 excess events 
corresponding to a significance of 7.1cr (Li & Ma 1983) are 
found at the nominal position of NGC 253. The overall statis- 
tics of MA with both sets of cuts are shown in Table Q] The 
best fit position of the source is RA 00 h 47 m 34.3 s +1.4 s , Dec 
-25°17'22.6"±0.3' (J2000), compatible at the < lcr level 
with the optical centre of NGC 253 at RA 00 h 47 m 33.1 s and 
Dec -25°17'18" (J2000). 

The squared angular distribution of y-ray candidate and 
background events relative to the position of NGC 253 as 
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TABLE 1 
VHE y-RAY statistics of NGC 253 



Cuts 


fp- 

°max 

(deg 2 ) 


N n 


N [T 


11 a 


Excess 


Significance 
(o-) 


MA (standard) 


0.01 


2240 


26224 


13.72 


329 ± 49 


7.1 


MA (faint) 


0.005 


571 


7816 


20.13 


183 ±24 


8.4 



Note. — Number of y- ray-like and background events along with the sig- 
nificance I Li & Ma 1983 1 and y-ray excess as obtained for the H.E.S.S. data 
using the MA method, a denotes the normalisation factor between signal and 
background exposure. 
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Fig. 2. — Distribution of ON events around NGC 253 and of OFF events 
from background control regions as obtained with MA, faint source cuts. The 
squared angular distribution has been produced using the Reflected back- 
ground model. Also shown is the point spread function of the instrument 
for this analysis as dotted line, assuming the y-ray emission originates from 
the optical centre of NGC 253. The normalisation of the model is adjusted to 
match the total y-ray excess in the range Odeg 2 - 0.01 deg 2 . 

shown in Fig. [2] is consistent with point-like emission. This 
constrains a potential source extension to less than 2.4' at 
3cr confidence level (cf. the extent of the starburst region 
in 12 CO(2-l) is about * 0.4' x 1.0'; |Sakamoto et al.pOlT 



In the standard picture, where the HE and VHE y-ray emis- 
sion from starburst galaxies originates from diffuse CR inter- 
actions, no variability of the y-ray signal is expected. The 
yearly light curve of the y-ray emission is stable over the four 
years of observations within errors. A fit of a constant flux to 
the yearly light curve yields a mean integral flux above 1 TeV 
of (8.5 ± 1.8) X 10" 14 cirT 2 s _1 with a^ 2 of 3.3 for 3 degrees 
of freedom and is stable over the four years of observations 

within errors. 

The H.E.S.S. data set previously published in |Acero et al.| 
(2009) comprises 119hrs of good quality data, collected in 
the years 2005, 2008 and 2009, and represents 2/3 of the data 
set used in this work. Based on this larger data set, a better 
determination of the spectral characteristics is now possible. 
All results presented in that publication are consistent with the 
findings presented here. 

2.5. Spectrum 

The differential energy spectrum derived from this data 
set is shown in Fig. [3] and Table [2] and is well described 
by a power law: dN/dE = Fq ■ (£/lTeV)~ r with photon 
index F = 2.14 + 0.18 sta t ± 0.30 sys and flux normalisation 
F = (9.6 + 1.5 sta , (+5.7, -2.9) sys ) x lO^TeV- 1 cirr 2 s" 1 
with a chance probability of 7%. This yields an integral flux 
above the energy threshold of 190 GeV of F(> 190 GeV) = 
(5.6 + 1.2) x 10~ 13 cm" 2 s _1 . With a flux of (0.21 + 0.05)% of 



the Crab nebula flux above the energy threshold, NGC 253 is 
the source with the lowest VHE y-ray flux detected so far. 

Since the VHE y-ray signal from NGC 253 is so weak, sys- 
tematic effects could potentially influence the spectral recon- 
struction. In-depth systematic checks on the background sub- 
traction and the spectrum calculation have been performed 
with the MA and BDT method. These tests have been used 
to estimate the systematic uncertainty of the VHE y-ray spec- 
trum presented here. For instance, the difference in the data 
selection procedures leads to a difference in live time of 12% 
and hence to a small difference in the data sets used to recon- 
struct the spectrum. Furthermore, the Earth's magnetic field 
bends charged particles in EASs and influences their devel- 
opment. This affects observables such as the shower image 
orientation in the Cherenkov camera, the image shape and 
hence the stereoscopic reconstruction. Although this effect is 
expected to be small, it is taken into account in the reconstruc- 
tion of the energy spectrum of NGC 253. In order to test the 
background systematics, harder selection cuts on y-ray like 
events have been performed, e.g. the faint cuts introduced in 
Section 2.3 and alternate background estimation techniques 
(such as the template background) have been studied as well. 
Precise comparison between the actual level of background in 
the whole field of view and predictions from the background 
model, excluding a circle of 0.25°around NGC 253, indicate 
that the background level is controlled at a level of * 1.5% 
over the full field of view, with the magnetic field introducing 
an azimutal asymmetry of ±2%. The magnetic field as well as 
potential systematic effects in the background determination 
mainly affect low-energy events. Therefore, additional tests 
have been performed, where only events with reconstructed 
energies above 0.6 TeV have been used in the spectral fit. 
All systematic tests have been performed with both analysis 
chains and resulted in a spread in normalisation between the 
two analyses of +60%, -30% and a variation in differential 
photon index of AF sys < 0.3. 

3. FERMI-LAT OBSERVATIONS AND DATA ANALYSIS 

Based on eleven months of data, NGC 253 has also been 
detected in the HE y -ray regime by the Fermi-LAT instrument 
(Abdo ~et al.| 2010a). It was later confir med as a HE y-ray 
emitter bas ed on 1 2 and 24 months of data ( |Abdo et al.|2 010b ; 



Nolan et al.|2012 1. The analysis and the results of a larger 30 
month data set are presented in the following. 

3.1. Fermi-LAT instrument 

The Fermi-LAT instrument is a pair-conversion telescope, 
capable of detecting y rays in the energy range between 
20MeV and 300 GeV. It consists of a tracker for the re- 
construction of the particle direction, a calorimeter which 
measures the energy of the incident particle, and an anti- 
coincidence system designed to suppress the charged-particle 
background. Data are normally recorded in survey mode, in 
which the whole sky is covered every two orbits. The instru- 
ment FoV is w 2.4 sr and it provides an angular resolution 
of < 1° at 1 GeV and < 0.2° at 10 GeV. A full description 
of the mission - and i nstrument-related details can be found in 
|AtwoodetaLlp009) . 

3.2. Data Set and Data Analysis 

The data set presented in the following comprises data 
from the commissioning of Fermi on 2008 August 4 th , (MJD 
54682) until 2011 February 3 rd (MJD 55595). The data anal- 
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Fig. 3. — Differential H.E.S.S. energy spectrum of NGC 253 as obtained with MA (shown as circles). Also shown is the LAT energy spectrum as obtained in the 
analysis of the 30 month s of data as de scribed in the text (shown as crosses). For both spectra lcr error bars are shown for spectral points and 95% upper limits 
according to Feldman & Cousins 1 199 8}. The s haded area represents the lcr confidence band from the simultaneous fit to the Fermi-LXt and H.E.S.S. data. Also 
shown are the predictions from Paglione et al. 1 1995) (solid grey), Domingo-Santamaria & Torres (2005 1 (dashed grey) and Rephaeli et al. 12010) (dash-dotted 
grey). 



TABLE 2 

H.E.S.S. SPECTRAL RESULTS OF NGC 253. 



Data 






Fo 


F(> £„,) r 




(GeV) 


(TeV 




( cm - s ) 


this analysis by MA 


190 


(9.6 ± 


1.5) x 10~ 14 


(5.6 ± 1.2) X 10 -13 2.14 ±0.18 


prev. analysis jAcero et al. 12009} 


220 






(5.5 ± 1.0) X 10 13 



Note. — VHE y-ray spectral results as shown in Fig.[3]for the MA. Th e phot on index T is derived from a fit of a power law to the spectrum. Only statistical 
errors are given. A comparison with the integral flux given in Acero et al. (2009) is also provided. 



ysis has been performed using events with reconstructed en- 
ergies between 200 Me V and 200 GeV and utilising the Fermi 
Science Tools (FST) package, version v9rl8p6rj Events 
of the "Diffuse" source class have been analysea using the 
P6JV3 -Diffuse instrument response functions (IRFs). Addi- 
tionally, events with zenith angles > 105° were excluded from 
the analysis due to a significant contribution of Earth-limb y 
rays^l 

The best fit position of NGC 253 has been determined by 
means of a maximum likelihood method, using the FST tool 
gtfindsrc. All events in a Region-of-Interest (Rol) of 10° 
around this best fit position have been used, and all sources 
within 15° were modeled to produce energy spectra and to 
calculate the test statistic (TS^jof the source. At the best fit 
position at RA 00 h 47 m 55.7 s and Dec -25°18'32.4" (J2000) 
(>95 = 5.7', 95% confidence level) a TS value of 105, corre- 
sponding to a statistical significance of a; 10 cr is found. 

Using the FST gtlike tool and an un-binned maximum like- 
lihood fitting procedure, the energy spectrum of NGC 253 has 

36 http: //fermi . gs£c .nasa.gov/ssc/data/analysis/ 

documentation/Cicerone/ 

J ' Note that as a cross check, the data have also been analysed with the 
more recent FST release, version v9r23pl and the P7Source_V6 IRFs, and 
found to give consistent results. 

38 Yts j s a measure of the significance and defined as TS = 
- 2 Alog( likelihood) between models including and excluding the source 
jMattox et al.|1996| . 



been derived. For this purpose, all sources listed in the Fermi- 
LAT 1-year catalogue in the 15° region around the best fit 
position were modeled with a power law {dN/dE cc E~ r ), 
with differential photon index T. In a second step, residual 
sources in the TS map in the Rol with a TS > 16 were in- 
cluded in the model as well. Compared to the 1-year and 2- 
year catalogue 20 and 16 more candidates, respectively, are 
included in the source model used in this work. In the min- 
imisation procedure, the integral flux in the energy range of 
interest and the photon index were left as free parameters for 
all sources in the Rol and fixed for all sources between 10 
and 15°. The diffuse Galactic and Extragalactic background 
components were modeled with the files gll_iem_vQ2.fit and 
isotropic JemjvQ2.tx^\ respectively. The energy spectrum is 
best described by a power law in energy with F = 2.24 + 
0.14 stat + 0.03 sys and an integral flux of F(0.2 - 200 GeV) = 
(4.9 ± 1.0 st at±0.3 svs )x 10~ 9 cirT 2 s _1 . The positional and spec- 
tral results of the likelihood analysis are summarised in Tab.[3] 
The flux in each energy band as summarised in Tableland 
shown in Fig. [3] has been reconstructed in the same way as 
for the full energy range, where the photon index and integral 
flux in that band were left as free parameters again. Note 
that all spectral points agree within 1 cr statistica l error with 
the spectral points and upper limits as reported in|Abdo et al. 



39 http : //fermi . gs£c . nasa . gov/ssc/data/access/lat/ 

BackgroundModels .html 
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TABLE 3 

Fermi-hAT spectrum and position of NGC 253. 



RA (J2000) 


Dec (J2000) 


r g5 


F(0.2-200GeV) T 


TS 


(deg) 


(deg) 


(deg) 


(HT'cm -2 s _1 ) 




11.982 


-25.309 


0.095 


4.9 ± 1.0 stat ± 0.3 sys 2.24 ± 0.14 stat ± 0.03 sys 


105 



Note. — Best fit position and results of the maximum likelihood analysis between 200 MeV and 200 GeV. 



TABLE 4 

i%rm(-LAT spectral results of NGC 253. 



TABLE 5 

Symbols, units and descriptions of quantities used in the discussion 



^Vnin ^max 


^NGC 253 


syst. uncer. 


(GeV) 


(10- 12 ergciri~ 2 s _1 ) 


(%) 


0.2 - 0.6 


1.9 ± 0.6 


10 


0.6-2.0 


1.7 ±0.3 


10 


2.0 - 6.0 


0.8 ±0.3 


15 


6.0 - 20.0 


1.2 ±0.5 


20 


20.0 - 60.0 


< 1.7 


20 


60.0 - 200.0 


< 8.7 


20 



Note. — HE y-ray spectral results of the maximum likelihood analysis in 
energy bands as shown in Fig. [3] lcr statistical errors and 95% upper limits 
are given. Syste matic uncertaintie s have been obtained using the bracketing 
method (see e.g. |Abdo et al.|2009) . 



(201C§I. 

The systematic uncertainty on the total spectru m has been 
estim ated following the bracketing method (e.g. Ab do et ah] 
2009), where the effective area is shifted up- and downward 
according to its systema tic uncertain t y in the corresponding 
energy range. Following Abdo et al. (2010b ), the systematic 
uncertainty on the flux as inferred by the systematic uncer- 
tainty of the effective area is of the order of 10% at 100 MeV, 
5% at 0.5 GeV and 20% at 10 GeV. As for the VHE y-ray light 
curve, there is no significant variability seen in the yearly HE 
y-ray emission. 

4. THE HE- VHE y-RAY SPECTRUM 

Within the lcr statistical errors, the H.E.S.S. measurement 
and the Fermi results are compatible, in respective value of 
photon index as well as respective normalisation. A simul- 
taneous, single power law fit to the H.E.S.S. and Fermi-LPtf 
spectral points results in a photon index of F s j m = 2.34 + 0.03 
and an energy flux at 1 GeV of (1.5 ±0.2) x 10~ 12 erg cirT 2 s" 1 
Fj This corresponds to an integral energy flux above 0.2 GeV 
of F™ eas (> 0.2 GeV) 5.3 x lO^ergcnrV. The fit has 

2 



&X of 8.27 for 7 degrees of freedom and a p-value of 30%. 
Note that the somewhat steeper index for the simultaneous 
fit compared to the individual measurements is fully compat- 
ible within the lcr statistical uncertainties of the Fermi and 
H.E.S.S. photon index. Even when taking into account a 30% 
downward shift to account for the systematic error in recon- 
structed normali satio n of the VHE y-ray spectrum as pre- 
sented in Section 2.5 the single power law fit has a p-value of 



The fit of a broken power law to the HE and VHE y-ray 
points with a break energy between 10 and 300 GeV (between 
H.E.S.S. and Fermi) results in a best-fit index change of AF = 
0.1 ± 0.3 with break energy at 300 GeV. Although, a photon 
index change between the Fermi and H.E.S.S. energy range 
can not be excluded, the hypothesis of no index change is 

40 Note that no forward-folding technique has been applied and the fit is 
performed only taking into account statistical errors. 



Symbol 


Units 


Description 


£ 


GeV 


energy of individual particle 


vm 
C PP 


GeV 


total non-thermal energy of pion- 




producing particles 


Q 


GeV" 'cm" 3 s _1 


differential volumetric particle source 
term 


Q 


GeVs-' 


volume-integrated particle source term 


Q" 


GeVs-' 


volume-integrated source term of pion- 
producing particles 



favoured given the available data. This result suggests that 
in order to explain the available data no spectral model which 
includes a break or turnover is needed. 

The integral (HE- VHE) flux above 200 MeV corresponds 
to a y-ray luminosity of l£(> 200 MeV) 7.8 x 10 39 ergs" 1 , 
for a distance of 3.5 Mpc. This is about a factor of ten larger 



than the y-ray luminosity L™ w (> 200 MeV) =* (7 - 10) x 



10 38 erg s 1 for the Milky Way galaxy, as estimated by Strong 
[etaLl ( [20T0T ). 

5. DISCUSSION 

The results presented in the previous sections have a num- 
ber of interesting implications for the nature and properties of 
the y-ray source in the starburst nucleus of NGC 253. Note 
that the discussion is based on the observationally favoured 
result that the spectrum of NGC 253 can be described by a sin- 
gle power law from the HE to the VHE y-ray regime. Quanti- 
ties used in the following are summarised in Table[5] This sec- 
tion is organised as follows: firstly the aspect of dominance of 
a hadronic scenario is discussed. Secondly, arguments for an 
energy-independent particle transport in the starburst region 
are given, the consequences of this picture are explored, and 
limits on the diffusion coefficient are presented. Thirdly, the 
y-ray flux estimated from the supernova rate is compared to 
the experimental results. Fourthly, the magnetic field in the 
starburst region is estimated and finally, the contribution of 
discrete sources to the y-ray luminosity of NGC 253 is as- 
sessed. 

5.1. Dominance of hadronic y-ray emission 

Early model predictions for the various contributions to the 
overall HE and VHE y-ray spectral energy distribution (e.g. 
Paglio ne et al. 1996|) and subsequ ent more detai l ed calcu- 
lations (e.g . Domingo-Santamarfa & Torres||2005[ [Rephaeli 
eTaLl|20T0l |Lacki et al.||2010| |2011| i, all consider diffuse y- 
ray emission and neglect the contributions of discrete sources. 
They are roughly compatible with the observational results. 
Three such model curves are shown in Fig. [3] 

These largely phenomenological models use primary 
electron-to-proton ratios similar to those inferred for the In- 
terstellar Medium (ISM) of the Milky Way galaxy and pa- 
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rameterise the diffusive escape times from the starburst re- 
gion, while assuming advection speeds of the order of sev- 
eral hundred kms" 1 . Also the secondary leptonic component 
from the decay of charged pions is included in the calcula- 
tion of the expected Bremsstrahlung and Inverse Compton 
y-ray emission. In detail the three model spectra in Fig. [3] 
are based on very different assumptions regarding the rela- 
tive magnitudes of the advective (r ac |) and energy-dependent 
diffusive (Tdiff) particle escape times on the one hand, and 
on the assumed source spectra on the other. These models 
assume Tdiff _ To(e p /(lGeV)~ os , where e p is the proton en- 
ergy. In Paglio ne et al.| ([T996) r a d > Tdiff for all energies ob- 
servecpj With an assumed proton differential source spectral 
index s = 2.2 the y-ray spectral energy density (SED) has a 
rather steep dependence oc 7 on the y-ray energy e y due 



to the dominance of diffusive escape. Rephaeli et al. (2010) 
obtain T a d > Tdiff. However, their source spectrum is assumed 
to be quite hard (s - 2.0). These authors model the emission 
from the galactic disk as a whole and find approximately a 
SED oc e~ 0,35 (see their Fig. 3). No spatial profiles are given. 
Domingo-Santamaria & Torres (2005) use for their main plot 
To = lOMyr and Tad = 0.3 Myr, resulting in T a d < Tdiff for 
£p < 1 TeV. Therefore, their SED should for low y-ray en- 
ergies be advection-dominated and correspond to the source 
spectrum. The corresponding curves in their Figs. 5 and 9 
indicate this source spectrum for an assumed particle source 
index s =s 2.3 up to e y equal to some hundreds of GeV. These 
curves indicate a significant softening beyond about 1 TeV. 
For a smaller value of to their SED should fall off strongly 
already beyond some considerably lower value 



-(s-2)-0.5 



Of 6y. 

The above model results show hadronically dominated y- 
ray emission up to several TeV, where the y rays are primarily 
produced in inelastic collisions between nuclear CRs and tar- 
get nuclei from the ambient ISM, and subsequent n° decay. 
Bremsstrahlung and Inverse Compton emission from primary 
and secondary electrons is not entirely negligible and in some 
of these models only by a factor of a few below the hadronic 
emission. This is in particular the case for y-ray energies be- 
low lOOMeV, where the 7r°-decay emission drops off, and at 
the high-energy end, where the harder Inverse Compton emis- 
sion in the Thomson limit can win over the 7r°-decay emission 
(if electrons at these energies do not experience significant ra- 
diative losses). Note however that the ratio of hadronic to lep- 
tonic emission depends on the assumed source spectrum, the 
form of Tdiff(e p ), and on the assumed electron-to-proton ratio 
(which has to be » than the canonical value of 1/100 that elec- 
trons dominate over the hadronic com ponent and overcome 
energetics problems, see e.g. Ohm & Hinton 2012). Therefore 
it is very likely that the hadronic emission dominates over the 
leptonic emission. That this is possible is basically a conse- 
quence of the very high gas density n„ in the starburst region. 
Aharonian et al. (2005 1 find a n„ — 580 cm -3 for a total gas 



mass o f 6 x 10' M , the value preferred by |Engelbracht et aL 
( 1998). Lower and higher gas masses as found by Mauers- 



berger et al. ( 1996[ ) and Sorai et al. ( 2000[ ), respectively, imply 
an uncertainty of this estimate of » 20%. Note that this esti- 
mate only takes into account the uncertainty in total gas mass 
and neglects any uncertainty in the starburst volume due to 
possible projection effects. Clearly this density is an average 

41 Due to a misprint in that paper the To in their Fig. 6 is in fact 1 Myr 
instead of the quoted tq = 10 Myr (T.A.D. Paglione, private communication). 



value - the actual gas density is probably extremely inhomo- 
geneous, given the large localised energy inputs from stars 
and in particular from their supernova explosions. Therefore, 
the local density for p-p interactions could differ from the av- 
erage density. However, it is unlikely that protons escape from 
the denser regions without significant losses. Adopting this n g 
corresponds to an average loss time due to inelastic p-p colli- 
sions of f pp « 1.1 x 10 s yr (this value is mildly dependent on 
the observed spectral index due to the energy-dependence of 
the p-p cross section <x pp ; a derivation of this value is given 
in the Appendix). Finally it has to be noted that f pp scales 
inversely linear with the inelasticity factor of the p-p collision 
(see Appendix). This can introduce an additional uncertainty 
on f pp of as 20%. 

We conclude that the y-ray emission from the starburst re- 
gion is likely to be dominated by hadronic interactions. In the 
following we will discuss the impact of system parameters 
such as SN rate, outflow velocity and particle diffusion co- 
efficient on the y-ray emission within the hadronic scenario. 
The possible additional rol e of discrete y-ray sources will be 
briefly discussed in Section 5.6 



5.2. Cosmic-ray escape 

In the hadronically dominated case, there are three main 
scenarios for the expected y-ray emission, depending on the 
(potentially energy-dependent) probability of cosmic rays es- 
caping the starburst region. In the case that the escape prob- 
ability is «: 1 at all energies, the system can be said to be 
calorimetric. The measured y-ray flux is significantly lower 
than the flux expected in the calorimetric limit for the canoni- 
cal parameters, su ch as the cosmic-ray acceleration efficiency 
as used in Section 5.4 Cosmic -ray escape must therefore be 



considered. There are two competing mechanisms to remove 
cosmic rays from the system: diffusion and advection. Domi- 
nance of advection (which is an energy-independent process) 
over the full energy range would result in a y-ray spectrum 
that approximately resembles the source spectrum. Diffu- 
sion on the other hand is an energy-dependent process and 
would lead to a spectral steepening of the source spectrum. 
At very low energies, the diffusion loss timescale Tdiff be- 
comes very long and eventually comparable to the advective 
loss timescale T a d. At the critical energy, where the diffusive 
transport takes over from advective transport and Tdiff = T a d, a 
spectral break in the source spectrum is expected. 

Empirically, for the Milky Way ga laxy, the measure ments 
of the CR secondary-to-primary ratio (Strong et al. 2007]) sug- 
gest a softening of the CR source spectrum by a factor oc E~ a , 
with 0.3 < a < 0.6, above a few GeV per nucleon as a result 
of energy-dependent diffusion. This is consistent with a com- 
paratively modest production rate of hot gas and CRs per unit 
area of the Galactic disk which allows most of the thermal gas 
to cool radiatively. At the same time the excitation of mag- 
netic field fluctuations by the escaping CR particles and their 
coupling to the thermal gas is equally modest. Although a 
CR-driven wind will develop at large distances from the disk, 
the observed particle spectru m in the Galact ic disk is domi- 
nated by diffusive escape ( jPtuskin et al.|1997] ). On the other 
hand, the wind from the starburst region in NGC 253 is very 
strong as a result of very strong gas heating that cannot be 
compensated by radiative cooling. The energy flux density in 
CRs driving magnetic fluctuations is expected to be very large 
as well, gi ven the high SN rate in the small star burst volume 
(see e.g. Aharonian et al.|2005| and Section [53] ). These non- 
linear effects massively diminish the role of particle diffusion 
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relative to global advection in the wind. Indeed, as shown in 
the next Section, the observed y-ray flux F™ eas (> 200 MeV) is 
well explained quantitatively by pure particle advection in the 
wind from the starburst nucleus, which req uires a CR produc- 
tion ©£sn - 10 5() erg per supernova (e.g. |Drury et al.|l989| 
for the starburst parameters of NGC 253. 

The diffusion time in NGC 253's starburst nucleus can be 
expressed as t^s = (H/2) 2 k~ 1 , where k denotes the diffu- 
sion coefficient and H m 60 pc is the height of the starburst 
region, positioned symmetrically to the galactic mid-plane. 
If fad 5 fdiff is required, then k should not exceed about 
3 x 10 27 cm 2 s _1 for all energies below the last H.E.S.S. flux 
point at 4.7 TeV that co rresponds to a CR energy of « 30 TeV 
(see e.g. Kelner et al. 2006 1. Such a small diffusion coef- 
ficient could be the result of strong wave excitation by the 
exceedingly c oncentrated CR pr oduction in the small star- 
burst volume (Aharonia n et al.||2005| l. This value can be 
compared to that of Bohm diffusion which should be consid- 
ered the slowest possible form of diffusion for a randomised 
magnetic field configuration. The Bohm limit is given by 
fBohm ~ 3 x 10 22 (e G ev/^G)cm 2 s _1 , where e Ge v is the CR 
energy in GeV and B^g is the magnetic field in fiG. For a 
particle energy of 30 TeV and a magnetic field of 100//G, as 
estimated in the next section, the Bohm diffusion coefficient 
is about 3 x 10 24 cm 2 s and thus still three orders of magni- 
tude smaller than the maximum diffusion coefficient deduced 
from the single power-law spectrum of the y-ray data. The 
upper limit to the diffusion coefficient estimated for the star- 
burst region in NGC 253 can also be compared to the average 
value inferred for the Milky Way gala xy which is given by 



* ga i « 1.5 x 10 3Q (E/1 TeV) 1 



(Ptuskin et al. 



1997 



Aharonian et al.|2005[ ). With such a comparison, extrapolat 



ing to a particle energy of 30 TeV, k < 5 x 10~ 5 A- ga i is found. 
This is an interesting requirement on the scattering strength 
of the magnetic field fluctuations^jin the starburst region, but 
certainly not an outrageous one, comparing with the Bohm 
scattering level. 

Dominant transport by advection is suggested by the obser- 
vation of a rather hard y-ray spectrum and the smooth align- 
ment of the spectrum in the HE and the VHE y-ray regime — 
the most noticeable observational result of this paper (see Sec- 
tion]?]). Therefore, taking the combined Fermi and H.E.S.S. 
spectra as single power law as indicated in Fig. [3] with differ- 
ential index r s ; m = 2.34, argues that energy-dependent dif- 
fusion is not important in the energy range covered rj In 
fact, advection alone already explains the magnitude of~F™ eas 
quantitatively. 

Therefore, two conclusions may be drawn here: (i) the 
observed spectrum is most likely the result of an energy- 
independent transport mechanism, i.e. advection and adia- 
batic expansion in the starburst wind and inelastic nuclear en- 
ergy losses, and (ii) the combined observed y-ray spectrum in 
the HE and VHE y-ray regimes might in this case correspond 
approximately to the mean spectrum of the ensemble of CR 
sources in the starburst region] 44 ] 

42 This value can also be compared to the gyroradius of a proton with 
1 TeV energy in a 100 fiG field, which is r g a 3.3 X 10 13 (£Tev/Bl00^G) cm at 
2A.U. 

43 Energy-independent diffusion can in principle occur as a result of large- 
scale turbulent motions in the gas and may also contribute to energetic- 
parti cle confinement in the gas that is systematically streaming with v w ind 
(e.g. Bykov 2001 Parizot et al. 20041. This effect is not included here explic- 
itly. 



5.3. y-ray flux estimate 

If, as argued in the previous section, CR escape from the 
starburst region is independent of particle energy, the ex- 
pected total y-ray energy flux F ™ p can be estimated simply 
from the parameters of the system. This independence is ap- 
proximately true also for the loss rate with respect to inelastic 
p-p collisions; in fact, this loss rate will be approximated by 
an average over the range of particle energies corresponding 
to the observed y-ray spectrum (see Appendix). In this case 
the particle transport equation can be simply integrated over 
the energy range of the pion-producing particles. I n a leaky 
box-type a pproximation for the starburst region (e.g. Berezin- 
skii et al]f l990) the result can in addition be integrated over 
the starburst volume. Given t hat the lifetime of the st arburst 
is about (2 - 3) x 10 7 yr (see |Engelbracht et al.|l998j > which 
is large compared to the advective loss time of about 10 5 yrs, 
the system is in a quasi-steady state. This results in the fol- 
lowing balance relation for the total non-thermal energy E pp 
of pion-producing particles in the starburst region: 



F" 



1 

Tad 



1 



1 



^"adiab 



(i) 



In the derivation of Eq.[T]the source spectrum and the result- 
ing spectrum of pion-producing particles in the starburst re- 
gion are both assumed to have approximately the form of a 
power-law spectrum in momentum oc p-( r ^ +2 \ as e xpected 
from t he theory of d iffusive shock acceleration (e.g. Bland- 
|ford & Eichle r|T987| see Appendix). 

In Eq. ]1J the quantity Q 71 — f„Q is the fraction f„ < 1, due 
to pion-producing particles, of the total input rate Q of non- 
thermal energy from the CR sources. Approximately f„ as 
3 - F s ; m = 0.66, assuming that the observed y-ray spectrum 
with F s j m = 2.34 resembles the source spectrum. The quantity 
Tad = (#/2)/v w ind ~ 10 5 y r denotes the advective los s time, 
where v w ; n d « 300kms _1 (Zirakashvili & Volk|2006| l is the 
velocity of the starburst wind at the top/bottom of the starburst 
region. The adiabatic loss time is given by r ac iiab in the accel- 
erating outflow. In a first approximation the flow speed V rises 
from zero at the galactic symmetry plane in the perpendicular 
direction z as \V\ = v w j n d[|z|/(i//2)] kms~', yielding an adia- 
batic loss rate l/r adiab * (r sim - l)VV/3 = (F sim - l)/(3r ad ). 
Finally, f pp (n g 0.5 c (crpp))" 1 , is an average energy loss time 
for inelastic, catastrophic proton-proton collisions, which is 
inversely proportional to the effective gas density n g . The 
inelasticity is taken as 0.5. With the mean cross section 
(°"pp) ~ 33 mb (see Appendix) this leads to a mean collisional 
energy loss time f pp 1.1 x 10 s (n g /580cm -3 ) -1 yr. 

The total hadronic y-ray energy flux density Fy Xp and E* p 
are connected by 



p exp 



And 2 '; 



(2) 



pp 



where 77 ^ 1/3 is the 7r°-fraction from overall pion-production 
through hadronic collisions. 

Eq. [2] assumes an optically thin y-ray emission region. In- 
deed the optical depth r yy for yy-absorption in the diffuse ra- 
diation field in the starburst region and the remaining part of 
NGC 253 is small compared to unity for the y-ray energies 



1 In fact, the -y-ray spectrum is slightly harder than the parent proton spec- tram I Kelner et al. 2006 Kamae et al. 2006 Karlsson 2008 1. 
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considered, i.e. t 77 < 0.1 for e y < 2TeV (see Inoue|20fT i. 



Assuming that the CR energy sources in the starburst region 
are the SN remnants, then Q = vsn®£sn, where vsn is the SN 
rate, < 1 is the CR production efficiency, and £sn is the to- 
tal hydrodynamic energy release per event. In order to comply 
with the overall energetics of CRs in the Milky Way galaxy, 
an average energy re lease Q£sn ~ 10 50 er g into nuclear CRs 
should be assumed (Dmr yet al.|1989| l, with an uncertainty of 
about a factor of 2. For £sn = 10 m erg this implies © « 0.1. 

It is of interest to compare this expected flux with the mea- 
sured flux F™ eas (> 200MeV) * 5.3 x 10~ 12 erg cnr 2 s _1 . For 

this purpose a nominal v alue of 0.03 yr 1 for v$n is assumed 
(Engelbracht et al. 19981. Inserting E* p from Eq. [Tjinto Eq. |2] 

leads tcE3 



Fy xp ^2.6x 10 



-SN 



11 -2 - 

erg cm s 



, 10 51 erg 





oTT 



VSN 



0.03yr- 



fn 

0.66 



2.6Mpc J 



(fpp/1.1 x 10 5 yr)/(r ad /10 5 yr) x l.ir sim + 3 



, 0) 



where all parameters like (r a d/10 5 yr) are written in terms of 
their nominal values. In the following the first bracket on the 
r. h. s. of Eq. [3]is referred to as the production term and the 
second bracket is referred to as the loss term of Eqj3] 

Considering nominal parameter values in Eq. BT the ex- 
pected flux Fy Xp 10" 11 ergcrrT 2 s _1 , and then F Y /F™™ » 
1.9. Given the uncertainties in the measurements of the nu- 
merous multi-wavelength parameters involved the two fluxes 
are quite close. This supports the general picture on which 
Eq.[3]is based. 

Trie question is, to which extent this result leads to physi- 
cally relevant bounds on physics quantities like the efficiency 
parameter ©£sn- There are two possibilities to reduce the ex- 
pected Fy Xp to the observed value F™ eas . The first of them can 
be achieved by reducing ©£sn compared to its value inferred 
for the Milky Way galaxy and/or by reducing the astrophys- 
ical parameter vs^/d 2 , i.e. by reducing the CR production 
term in Eq. [3] The second possibility is to increase the ratio 
Tpp/T a d either by decreasing the effective gas density n g seen 
by the pion-producing particles during their advective escape, 
relative to the average gas density observed, or by increas- 
ing v w ind- In this second case the source term could be kept 
at its nominal value, especially 0£sn could be kept at the 
value of 10 50 erg. Decreasing f pp /T a d relative to its nominal 
value, on the other hand, is quite implausible for two reasons: 
firstly, the wind velocity cannot become significantly smaller 
without causing difficulties to explain the spat ial extent of the 
starburst region as seen in radio observations (Zirakashvili & 
|Volk|[2006| . Secondly, the total gas density cannot signif- 
icantly increase, given the observations of the molecular gas 
mentioned earlier. This means that the loss term in Eq.[3]has a 
maximum value for the default parameters in Eq. [3] As a con- 
sequence the product vsn©£sn should at most decrease by a 
factor 1/1.9, if at all, from its value 10 50 erg x 0.03 yr -1 for the 
nominal parameters. Effectively, 0£sn depends on the value 
vsn- Only for vsn » 0.03 yr -1 should 0£sn on average be 



45 The assumption that the parent proton spectrum and the resulting y-ray 
spectrum have the same power-law index infers an error of 
(see also Appendix). 



20% in E" pp 



small compared to 10 50 erg, i.e. in contrast to the average sit- 
uation in the Milky Way galaxy and theoretical expectations 
for individual supernova remnants there. 

To derive the quantities of NGC 253 the reference distance 
of d — 2.6 Mpc has been used most widely in the literature. 
However, as discussed in the Introduction, this distance has 
recently been revised to d — 3.5 Mpc. Therefore, an appro- 
priate scaling of the astronomical parameters in Eq. [3] needs 
to be considered. Indeed, the total gas mass as determined 
from the CO line flux scales as d 2 , whereas the starburst vol- 
ume scales as d 3 . This implies that the gas density scales as 
d~ l and f pp cc d. The supernova rate is derived from the FIR 

continuum flux and scales as vsn 00 d 2 . The wind velocity is 
derived to be consistent with the geometry of the radio bright- 
ness distribution yielding v w j„d cc d. As a consequence r a d is 
independent of d. For d — 3.5 Mpc the nominal value of F y xp 
in Eq. [3] is therefore about 8.3 x 10~ 12 ergcm -2 s _1 . Follow- 
ing the above physical arguments and changing the distance 
to d = 3.5 Mpc, the production term on the r.h.s. of Eq.[3]must 
not become smaller by more than a nominal factor 1 /T.6, in 
order to reduce F 7 xp to F ™ eas . 

In general, the simple model presented here agrees quite 
well with the observed values, given the observational uncer- 
tainties of the astronomical parameters and the possibility that 
for SN remnants in such a dense medium the total non-thermal 
energy generated per event 0£sn might indeed be lower by a 
factor a 1/1.6 than typically assumed for an object in the av- 
erage ISM of the Milky Way galaxy. 

5.4. Hadronic calorimetry 

In this scenario, to find the extent to which the starburst re- 
gion behaves calorimetrically in the presence of advective and 
diffusive escape, the total energy production rate in hadronic 
collisions L co u = And 2 F™^ I r] is compared with the total pro- 
duction Lcr{tt) = /,tVsn0£sn of CRs capable of producing 
hadronic y rays. Suc h a comparison for NGC 253 has previ- 
ously b een made by|Aharonian et al.| ( [2005| l; |Loeb & Wa xman 
(2006), and Thomps on et al.| ( |2007| l. LcbSjt) depends on the 
fraction of energy available for pion-production f n a 3 - s, 
with s being the source spectral index. In the present scenario 
fn - r s ; m and therefore f„ = 0.66. 

Ado pting a value of 0.03 yr -1 for vsn (Engelbracht et al. 
1998| l and assuming for ©£sn ~ 10 5() erg again results in 



^coll 

LcrM 
Esn 



0.21 x 



rmeas 
7 



5.3 x 10~ 12 ergcirr 2 s" 



10 5I erg 



© 
01 



0.03yr- 



0.66 



2.6Mpc 



(4) 



This fraction of about 20 percent is then a measure of the ex- 
tent to which the starburst region is calorimetric with respect 
to its hadronic interactions, dissipating its own non-thermal 
output; and since vsn is proportional to the FIR luminosity 
and hence to d 2 , this result is formally independent of the dis- 
tance d. The ratio L r f ,u / 'Lrp (n) is comparable to the value 



found by 



Lacki et al. 



"(201 l|l. 



Like Fy v /F^\ it is affected 
by the uncertainties of tEelnput parameters. Adopting a larger 
ratio of v^/d 2 would decrease this value, whereas a lower ef- 
ficiency ©£sn would increase it. The difference between F™ p 
and F™ eas , discussed above, shows that the value of the pro- 
duction term of Eq. [3] should possibly be smaller than unity 
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but not smaller than a 1/1.6 for d — 3.5 Mpc. This suggests 
that the calorimetric fraction may be larger than 20 percent, 
possibly reaching up to about 30 percent. 

5.5. Cosmic-ray energy density and magnetic field strength 

Within this framework the non-thermal energy density {/* 

of the 7T° -producing particles in the starburst region can be 
simply calculated from Eq.|2]by substituting Fy Xp = ,F™ eas and 
dividing the resulting E* p by the estimated volume. This gives 

U pp w 230eVcirT 3 , independent of d. This is a lower limit 
for the total non-thermal energy density U pp . Since it is diffi- 
cult to estimate the contribution of the lower-energy particles 
produced in the sources, because of their poorly determined 
ionisation losses in the outflow, one can only give an upper 
limit Upp < Upplf„ ~ 340eVcirT 3 . The values in this range 
are more than a hundred times larger than the CR energy den- 
sity in the Milky Way galaxy and, by implication, on average 
in the disk of NGC 253 (see Introduction). Assuming, as for 
the Milky Way, a CR scale length of about 1 kpc for the ex- 
tended disk of NGC 253 in comparison to the w 40 pc gradient 
scale in the starburst region (for d — 3.5 Mpc), the latter's CR 
pressure gradient is almost 4 orders of magnitude larger than 
that on average in the disk. 

Equating U pp with the magnetic field energy density B 2 /87r 
in the sense of a conventional equipartition argument for mag- 
netic field and CRs, results in a range for the r.m.s. magnetic 
field strength 90juG < B eq < 120 fxG. Note that the estimates 
of Upp and B eq are also affected by the uncertainty of the in- 
elasticity factor for p-p c ollisions that could be w 20% lower 
(Musulmanbekov] |2004| l. This is to be compared with the 
(160 ± 20) fiG equipartition field st rength for the so-calle d nu- 
clear region, estimated recently by Heesen et al. (2011 1 from 
radio continuum observations and assuming a proton to elec- 
tron ratio of 100. 

5.6. Contribution of discrete sources 

It has recently been suggested that the HE y-ray emission 
from NGC 253 might indeed come from interactions of dif- 
fuse nuclear CRs with the ambient gas, but that this contribu- 
tion should substantially diminish at energies > lOGeV as a 
result o f energy-dependent diff usion losses from the starburst 
region (Mannheim e t al.|2012) . The observed TeV emission 
was argued to be due to Inverse Compton emission from unre- 
solved pulsar wind nebulae (PWN) instead, i.e. from a sepa- 
rate population of discrete sources. For such PWNe the evolu- 
tion was assumed to be the same in the disk of the Milky Way 
galaxy and in the extremely high-density environment of the 
starburst nucleus of NGC 253. It has however to be noted that 
electron cooling times in typical starburst environments are 
very short and electr ons are expected to c ool on timescales of 
a few hundred years (Ohm & Hinton 20T2}. 

While in the present paper it is argued that advective re- 
moval should dominate over diffusive losses of energetic nu- 
clear particles deep into the TeV range, a discussion of the 
role of PWNe in such an environment is beyond the scope of 
the present work. On the other hand, the observations show 
that the spectrum of NGC 253 can be well described by a sin- 
gle power law over the combined HE and VHE energy range, 
disfavouring two distinct spectral components and suggesting 
that the same physical processes dominate over the entire en- 
ergy range in question. Therefore, additional discrete y-ray 
sources seem to play a minor role. 



6. SUMMARY AND OUTLOOK 

The analysis results of 177 hours of H.E.S.S. data obtained 
in observations of the starburst galaxy NGC 253 are reported. 
The reconstructed energy spectrum is best described by a po 
wer law with differential photon index F = 2.14 ± 0.18 sta t ± 
0.30 sys and differential flux normalisation at 1 TeV of Fq = 
(9.6 ± 1.5 sta t (+5.7, -2.9) sys ) x 10~ 14 TeV" 1 cnr 2 s _1 . In addi- 
tion to the H.E.S.S. data, the analysis of the 30 months Fermi- 
LAT data set revealed an improved best fit position compati- 
ble with the optical centre of the galaxy and with the H.E.S.S. 
source within statistical errors. The reconstructed differential 
photon index is T — 2.24 ± 0.14 stat + 0.03 sys and the integral 
flux between (0.2 - 200) GeV is >(0.2 - 200 GeV) = (4.9 ± 
1.0 st a t ±0.3 sys ) x 10~ 9 cnr 2 s _1 . The HE and VHE y-ray spec- 
tra of the starburst region can be described by a simultaneous 
power-law fit with differential photon index r s j m = 2.34±0.03 
and a fit probability of 30%. This result implies that no spec- 
tral break or turnover is required to explain the y-ray data. 
The corresponding total energy flux density corresponds to 
F E {> 200 MeV) * 5.3 x 10~ 12 erg cm" 2 s _1 . Assuming the 
remaining disk of NGC 253 to be quantitatively similar to the 
Milky Way galaxy, the starburst region outshines the rest of 
NGC 253 by an order of magnitude in HE and VHE y rays, 
consistent with the detection of the object as a H.E.S.S. point 
source. 

Model predictions which assume a dominantly hadronic 
origin of the y-ray emission are roughly compatible with the 
spectral results presented in this work. For a set of reason- 
able parameters the CR energy, which is lost in p-p interac- 
tions and partly re-appears in 7r°-decay y-ray production, is in- 
ferred in the present work as w 20% and possibly up to as 30% 
of the total non-thermal energy produced in the starburst re- 
gion, assuming a distance of 3.5 Mpc and a 10% efficiency 
for cosmic -ray acceleration in starburst supernova remnants. 
Note however, that the multi-wavelength observables are only 
known within a considerable error margin, which can change 
these percentages significantly. CRs are also removed by dif- 
fusion and advection from the starburst region. Since the for- 
mer process is energy dependent, a spectral steepening with 
energy would be expected. The smooth alignment of the HE 
and VHE y-ray spectra over four decades in energy hence in- 
dicates that advective losses in NGC 253 most likely dominate 
from a few GeV to more than 10 TeV. Even at such high en- 
ergies the diffusion coefficient would still be more than two 
orders of magnitude larger than the Bohm diffusion coeffi- 
cient. It therefore seems likely that the observed spectrum 
can be characterised by the same photon index as the average 
particle accelerator in the starburst region. 

The form of the y-ray spectrum of NGC 253 can be com- 
pared with another starburst galaxy, M 82 in the Northern 
Hemisphere, det ected by the Fermi-L AT ( Abdo et al.|2010a i 
and VERITAS ( jAcciari et al.|[2009] > collaborations, respec- 
tively. If one looks at the corresponding HE and VHE y-ray 
spectra of M 82, the overall shape looks rather similar to the 
one presented in this paper for NGC 253. Even though the 
starburst in M 82 is in all probability triggered by the interac- 
tion with the companion galaxy M81, the spectral similarity 
is consistent with the assumption that the CR sources in both 
galaxies produce similar energetic particle spectra and may 
therefore be of the same nature. 
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APPENDIX 



Eq. [T|of the main text can be derived from the transport equation for the isotropic part f(x*,p,t) of the particle momentum 
distribution in a simple model, negle cting the speed of the scatte ring fluctuations compared to the fluid mass velocity V as well 
as diffusion in momentum space (e.g. Blandford & Eichler| 1 987) 1 : 



at p op V 3 / r nn 



/ 



(Al) 



\_d_ 

p 2 dp\ ^ , . P p 

where k(j>) denotes the spatial diffusion coefficient, f pp is the average energy loss time as a result of inelastic, catastrophic nuclear 
collisions (assumed to be energy-independent for p > 2m p c and infinite for p < 2m p c, see below), and Q(j?,p,t) is the particle 

production rate of the CR sources. The adiabatic loss rate of a particle in the accelerating flow is given by p/3VV. In the 
energy range of pion-producing nuclei ionisation losses are neglected. Assuming a steady state {df/dt = 0) and neglecting 
particle diffusion (k = 0), the flow velocity V is approximated to be perpendicular to the galactic disk mid-plane, varying as 
V = (0, 0, v w i n d X z(H/2)~ l ) for < \z\ < H/2, where the constant parameter v w i nd is the wind velocity at the starburst boundary. 
To obtain the balance relation for the kinetic energy of the pion-producing CRs E pp — J y (fx An Jj" 1 ™ dpp 2 E\^ n f in the 



starburst region, Eq. Al is multiplied by the particle kinetic energy E^ n and integrated ov er the relevant mom enta as well as over 



the spatial starburst volume Vsb ■ m a leaky box approximation for the starburst region (e.g. Berezinskii et al.|1990 1, f(x, p) and 



' pp 



are assumed to be spatially uniform. It is then clear for this particle distribution that it has the same momentum dependence as the 
spatially averaged source production rate f. d 3 xQ(x, p)/Vsb for » m j n < p < p max . In addition, it is assumed that the production 
spectrum in the sources is not strongly influenced by energy-dependent losses in the sources themselves. The averaged particle 
source production rate can then be assumed to have a power-law dependence oc p-( s+T > in particle momentum p for all momenta 



above the i njection energy, consistent with an origin of the CR particles from diffusive shock acceleration (e.g. Blandford & 
Eichler|1987]l. And s as F s j m , where r s j m = 2.34 is the power-law index of the observed y-ray flux. The minimum momentum 
roughly corresponds to GeV energies. Here p 



in Eq. 



Al 



is taken as 2m p c. Since the maximum observed y-ray energy is 
4.7 Te V, corresponding to a proton energy of cp max as 30 TeV, momenta p > /? max make a negligible contribution. Therefore, it is 
possible to put p m . dx = oo. 

In order to obtain an analytical estimate the above integral for E pp is approximated by using the relativistic formula E^ n = cp, 
for p > 2m p c ( |Drury et al.[T989] ). 
Integrating the individual terms in Eq. |Al| by parts over momentum results in : 



F" 
£ pp 



1 



1 



1 



where r ad = (H/2)/v wind , 

^"adiab 



^"adiab ^~ 

3r ad /(r sim - 1), and Q" * 4n £ 



(A2) 



2m c 



dpp 2 Q. This quantity Q K - f„Q defines the fraction 



f n < 1 of the total particle energy production Q = J d 3 x4n C° dpp 2 E^ n Q, integrated over all momenta, that is available for 
pion production. Consistent with the approximation for Q 71 the integral for Q is approximated by using the newtonian formula 



£kin = p /2m p for p < 2m v c and the relativistic formula E^ = cp for p > 2m p c ( [Drury et al 



1989). 



Assuming 2 < s < 3, p m ^ <s 2m p c, for the injection momentum p m ^, and /? max » Zm v c, then gives: 



J2m v c 



dpp 2 cp(p/2m p c) 



-s-2 



4ttQ c/(s - 2)(2m p c) 4 



and 



Q - Q" * 4ttG 



Pom 



dpp 2 p 2 /2m p (p /2m p c) 



■s-2 



4nQ c/(3 - s)(2m p c) 4 



(A3) 



(A4) 



Therefore, Q" at (3 - s)Q and /j*3-j for the assumed m omen tum spectrum. For s - r s j m = 2.34, f„ as 0.66. For a very hard 
source spectrum with s close to 2, f n would be close to 1. Eq. |A2| corresponds to Eq.[T]of the main text. 
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The approximate quantity f pp , i.e. the mean loss time that appears in Eq. Al arises from averaging the energy loss rate 
r pp 1 (£') over the differential energy spectrum of the colliding particles. Approximating the differential energy dependence of 
the particle spectrum by that of the observed y-ray spectrum implies that the average can be taken over a power law spectrum 
/(^kiii) 00 E kin s,m which has the same index r s ; m = 2.34 as the differential y-ray flux as a function of the y-ray energy e r in the 
energy region observed for NGC 253, thus defining: 

J_ = fdEEf(E)/T pp 
r PP J dEEf(E) 

The integration limits are the particle (proton) energies corresponding to the observed energy range of the y-ray energy spectrum. 
Since T p l p — n g c/i n cr pp (£ , kin ) with the total cross section cr pp for inelastic pp-collisions and inelasticity factor f m ss 0.5 , cf. 



Aharonian ( 2004| l, this amounts to the following average over the inelastic nuclear cross-section <x pp : 

1 "gc/in JdEE <t w (E) f(E) 



T 



l>l> 



fdEEf(E) 



(A6) 



It is a convenient approximation to choose the integration limits as EJ™ = 1 GeV, and as £J™ X ~ (4.7 TeV/0.17) * 30TeV 
according to the maximu m y-ray energy e y = 4.7 TeV observed. For E^ n > 1 GeV the cross section <x pp is approximated in the 
following form (Aharonian 2004): 



(E m ) ~ 30 [0.95 + 0.06 ln^/GeY)] mb (A7) 



The corresponding integrals can be calculated analytically. Since the spectrum is relatively soft, the resulting effective cross 
section (<x pp ) as 33 mb corresponds to a rather low effective particle energy (E^ n ) as 12 GeV. Due to the slightly softer source 
spectral index compared to the resulting y-ray spectral index (s - 2.4 vs. r s ; m = 2.34), the assumption of equal spectral indices 
infers an error of « 20% of the proton energy density E" pp in Equations |3] and |4] of the main text (S.R. Kelner, priv. comm.). This 
is certainly a negligible effect, compared to all other uncertainties. 
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